/**
  SUNNY Project, Anyang Normal University, IMP-CAS
	\class TAMBSDList
	\brief A list to store many-body Slater determinants, classified by M, which is
	the 3rd component jz of the total angular momentum. So M is the same for each
	member of this list.
	\author SUN Yazhou
	\since 2020/01/31
	\date Last modified: 2022/12/05 by SUN Yazhou
	\copyright 2020-2023 SUN Yazhou
	\copyright SUNNY project, Anyang Normal University, IMP-CAS
*/

#include <iostream>
#include <algorithm>
#include <iomanip>
#include "config.h" // contains user switch of definition of ROOT and OpenCV
#include "TAMBSDList.h"
#include "TAMBSD.h"
#include "TAException.h"
#include "TABit.h"
#include "TASPStateManager.h"
#include "TAMBSDManager.h"

using std::cout;
using std::endl;
using std::setw;

class TAMBSD;

TAMBSDList::TAMBSDList(){
#ifdef _ROOT_
  fMBSDTree = new TAMBSDTree("mbsdList", "Many-Body Slater Determinant Tree");
#endif
} // end ctor

TAMBSDList::~TAMBSDList(){
#ifdef _ROOT_
  if(fMBSDTree) delete fMBSDTree;
#else
  for(auto &p : fMBSDVec) if(p){ delete p; p = nullptr; }
#endif
} // end the destructor

bool TAMBSDList::MBSDFilter(ull bit) const{
  static auto bm = TAMBSDManager::Instance();
  return bm->MBSDFilter(bit);
} // end member function MBSDFilter

void TAMBSDList::Add(ull bit){
#ifdef _ROOT_
  if(MBSDFilter(bit)) fMBSDTree->Fill(bit);
#else
  if(MBSDFilter(bit)) fMBSDVec.push_back(new TAMBSD(bit));
#endif
} // end member function Add

void TAMBSDList::Add(int *spsArr, int nParticle){
  TABit bit; bit.SetBit(spsArr, nParticle);
  Add(bit.to_ulong());
} // end member function Add

ull TAMBSDList::Bit(ull i){
  if(i < 0 || i > GetNBasis())
    TAException::Error("TAMBSDList", "Bit: mbsdIndex %d out of range.", i);
#ifdef _ROOT_
  return fMBSDTree->GetMBSDInBit(i);
#else
  return fMBSDVec[i]->Bit().to_ulong();
#endif
} // end member function Bit

void TAMBSDList::Print(){
  int n = GetNBasis();
  if(n > 500) n = 500;

  cout << setw(5) << "index" << setw(3) << "pi";
  cout << setw(3) << "2M" << setw(3) << "2T" << setw(5);
  cout << "E" << setw(8) << "bit" << setw(16) << "ManyBodySD Arr" << endl;
  for(int i = 0; i < n; i++){
    cout << setw(5) << i;
    TAMBSD((*this)[i]).Print();
  } // end for over i
} // end member function Print

/// Print all the mbsd-s in bit mode
void TAMBSDList::PrintInBit(){
  int n = GetNBasis();
  if(n > 500) n = 500;
  cout << "ManyBodySD Array in bit" << endl;
  for(int i = 0; i < n; i++) TAMBSD((*this)[i]).PrintInBit();
} // end member function PrintInBit

ull TAMBSDList::GetNBasis() const{
#ifdef _ROOT_
  return fMBSDTree->GetEntries();
#else
  return fMBSDVec.size();
#endif
} // end member function GetNBasis

// note that index is the mbsd index; i is the i-th particle
TASPState *TAMBSDList::GetSPState(ull mbsdIndex, int i){
  static vector<TASPState *> &spv =
    TASPStateManager::Instance()->GetSPStateVec();
  if(i >= int(spv.size()))
    TAException::Error("TAMBSDList",
      "GetSPState: Subscript of the requested sp orbit is out of range.");

  TABit b(Bit(mbsdIndex));
  int np = -1, j, n = b.size();
  for(j = 0; j < n; j++){
    if(b.test(j)) np++;
    if(np == i) break;
  }
  TASPState *sp = spv[j];
  if(!sp) TAException::Error("TAMBSDList", "GetSPState: Required sp is null.");
  return sp;
} // end member function GetSPState

/// \retval <rr|a+_p * a_q|cc>
int TAMBSDList::Integral(int rr, int p, int q, int cc){
  TABit rBit((*this)[rr]);
  rBit.Annihilate(p); if(!rBit.GetPhase()) return 0;

  TABit cBit((*this)[cc]);
  cBit.Annihilate(q); if(!cBit.GetPhase()) return 0;

  return rBit*cBit;
} // end member function Integral(rr,p,q,cc);

/// \retval <rr|a+_p*a+_q * a_s*a_r|cc>
int TAMBSDList::Integral(int rr, int p, int q, int s, int r, int cc){
  TABit rBit((*this)[rr]);
  rBit.Annihilate(p); if(!rBit.GetPhase()) return 0;
  rBit.Annihilate(q); if(!rBit.GetPhase()) return 0;

  TABit cBit((*this)[cc]);
  cBit.Annihilate(r); if(!cBit.GetPhase()) return 0;
  cBit.Annihilate(s); if(!cBit.GetPhase()) return 0;

  return rBit*cBit;
} // end member function Integral(rr,p,q,r,s,cc);
/// \retval <rr|a+_p*a+_q*a+_r * a_u*a_t*a_s|cc>
int TAMBSDList::Integral(int rr, int p, int q, int r, int u, int t, int s, int cc){
  TABit rBit((*this)[rr]);
  rBit.Annihilate(p); if(!rBit.GetPhase()) return 0;
  rBit.Annihilate(q); if(!rBit.GetPhase()) return 0;
  rBit.Annihilate(r); if(!rBit.GetPhase()) return 0;

  TABit cBit((*this)[cc]);
  cBit.Annihilate(s); if(!cBit.GetPhase()) return 0;
  cBit.Annihilate(t); if(!cBit.GetPhase()) return 0;
  cBit.Annihilate(u); if(!cBit.GetPhase()) return 0;

  return rBit*cBit;
} // end member function Integral(rr,p,q,r,s,t,u,cc);
